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We describe analytically the cosmological evolution of primordial gravity waves (tensor pertur- 
bations) taking into account their coupling to relic neutrinos. We prove that as a consequence of 
causality the neutrino-induced phase shift of subhorizon tensor oscillations tends on small scales to 
zero. For the tensor modes that reenter the horizon in the radiation era after neutrino decoupling 
we calculate the neutrino suppression factor as 1 — 5p„/9p + 0[(p„/p) 2 ]. This result is consistent 
with the value obtained for three neutrino flavors by Weinberg and is in agreement with numerical 
Boltzmann evolution. A minimal formula with the correct asymptotic form on small and large scales 
reproduces to about ten percent the evolution on all scales probed by the CMB. A more accurate 
solution (in terms of elementary functions) shows that the modes reentering the horizon in the 
radiation era are slightly enhanced and the phase of their temporal oscillations is shifted by sub- 
dominant nonrelativistic matter. The phase shift grows logarithmically on subhorizon scales until 
radiation-matter equality; the accumulated shift scales for k S> fc oq as lnfc/fc. The modes reentering 
the horizon after equality are, in turn, affected by the residual radiation density. These modes follow 
the naive matter-era evolution which is advanced by a redshift and scale independent increment of 
conformal time. In an appendix, we introduce a general relativistic measure of radiation intensity 
that in any gauge obeys a simple transport equation and for decoupled particles is conserved on 
superhorizon scales for arbitrary initial conditions in the full nonlinear theory. 



I. INTRODUCTION AND OVERVIEW OF 
RESULTS 

The primordial power of gravity waves, or tensor met- 
ric perturbations, is a direct probe of the inflation energy 
scale [1-4], unknown at present. Observation of their B- 
mode signature [5, 6] in polarization of the cosmic mi- 
crowave background (CMB) is a target of several devel- 
oped experiments. Gravity wave modes are gcncrically 
frozen while the Hubble scale is smaller than the mode 
wavelengths. Yet, a mode starts to evolve as soon as it 
"enters" the Hubble horizon. This evolution determines 
the pattern of the tensor signatures in the CMB spec- 
tra. The modes processed by the evolution are to be 
constrained by ground and space-based interferometers, 
resonant mass detectors, or timing of millisecond pulsars. 

While numerical calculations are helpful in theoretical 
studies and may be necessary for data analysis, exact an- 
alytical results provide a framework for our understand- 
ing of physical phenomena. Quantitative analytic de- 
scription of gravity wave evolution in the realistic cosmo- 
logical setting, roughly established today, is the subject 
of this paper. A consistent, involving no fitting and ad 
hoc approximations, study in elementary formulas with 
clear interpretation can be performed to several percent 
accuracy. 

We invoke both the Fourier and real space views of cos- 
mological dynamics. The advantage of Fourier decompo- 
sition is decoupling of modes with different wavenumbers 
in linear order. There are also merits in the complimen- 
tary real space approach. First, the real space formula- 
tion manifests the locality of interactions and causality 
of solutions. Causality imposes severe restrictions on the 



admissible evolution, and real space is a preferred start- 
ing ground to explore the implications. Second, only in 
real space we have an appealing formulation of the full, 
non-linear, cosmological dynamics. Understanding of the 
linear regime in real space is thus essential for connecting 
linear and non-linear theories. 

Evolution of a tensor (transverse traceless, helicity ±2) 
mode hij of the metric (12) with a wavevector k obeys a 
wave equation, e.g. [7-9], 

hij + 2Uhij + k 2 h l3 = 167rGa 2 Sy. (1) 

The mode dynamics is affected by the Hubble expan- 
sion rate H. = a/a and the tensor component of mat- 
ter anisotropic stress £,j (20). After the horizon entry 
the mode undergoes continuous oscillations, which are 
damped by the Hubble expansion. On subhorizon scales 
{k 3> H), the damping becomes adiabatic, h cx 1/a, and 
the oscillation period settles at 2ir/k. The asymptotic 
amplitude and phase of the subhorizon oscillations de- 
pend on the mode evolution during the horizon entry. 

The modes with k > fc cq = r" 1 ~ 10 _2 Mpc _1 en- 
ter the horizon when photons and neutrinos make up a 
substantial fraction of the overall energy density. This 
applies to all the modes accessible to direct detection 
and to the modes exciting the CMB multipoles with 
I ^ ^ydecAeq ~ 130. These modes also contribute appre- 
ciably to the lower multipoles of CMB polarization up to 
the rise of the reionization bump at I ~ 10. 

If anysotropic stress were negligible then during radia- 
tion domination /i( rad >°) = jo((p) = simp/ (p. Here and for 
the entire paper 

ip = kr (2) 
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and mode normalization is h — > 1 as a — > 0. In real 
space, the neutrinoless radiation era evolution of grav- 
ity waves is described by a top-hat Green's function 

/l(rad,0) = J-6{T-\x\). 

However, recently Weinberg [10] showed that the con- 
tribution of free-streaming neutrinos to anisotropic stress 
and the consequent gravity wave damping [8, 11] is by 
no means tiny. The subhorizon amplitude of tensors and 
the induced by them CMB polarization were found sup- 
pressed by a factor A w 0.8 [10] on all scales entering 
the horizon from neutrino decoupling to matter domina- 
tion, H~\ cc ~ 0.1 kpc to k~^. All tensor power spectra 
are suppressed on these scales by A 2 ,. 1 

Weinberg's result was derived by calculating the 
anisotropic stress E of free streaming neutrinos as a line 
of sight integral over the metric perturbation h. Hence, 
cq. (1) was converted into a closed integro-differential 
equation. The source on its right hand side, however, 
was a complicated integral, which depended on hij(r) 
nonlocally and should be fully computed numerically at 
every evolution step. The amplitude and phase of the 
subhorizon oscillations could then be found as the ulti- 
mate output of numerical evolution with this source from 
superhorizon to subhorizon scales. 

A real space approach provides an instant proof that 
regardless of the abundance of neutrinos, or presence of 
additional exotic species, the phase of subhorizon tensor 
oscillations in the radiation era is unshifted with respect 
to the zero-stress solution ft,( rad '°): 

h(rad) _> A smip/ip (3) 

as ip — > oo. 2 By causality of tensor metric perturbations 3 , 
their Green's function h(r, x) which is localized at the 
origin at r — > should remain zero beyond the particle 
horizon: 

h(r, x)=0 for \x\ > t. (4) 



1 Photon anisotropic stress, which develops after the recombina- 
tion, additionally suppresses a narrow range of tensor modes with 
k ~ 0.005 Mpc- 1 by up to 7%. The affected arc the modes that 
enter the horizon after photons decouple but before the dynam- 
ical role of radiation becomes entirely negligible. This effect, 
however, has little consequence for CMB polarization, generated 
earlier during the decoupling, and even for the tensor ISW im- 
pact on the CMB temperature: C ; ss and the tensor component 
of Cf T are reduced by less than 0.5% and than 1% respectively. 

2 An analogous result exists in the scalar sector [12]. On small 
scales, the phase of photon-baryon acoustic oscillations can be 
shifted only by gravitational coupling to species whose pertur- 
bations propagate faster than the photon-baryon sound speed. 
Among such species are the free-streaming neutrinos and a hypo- 
thetical scalar field with non-negligible density during radiation 
era (early quintessence). 

3 Propagation of gauge-dependent scalar and vector perturbations 
of the metric, density, or a field generally is not causal, even if 
there is no physical violation of causality. Nevertheless, when 
the initial conditions are adiabatic, scalar perturbations expand 
manifestly causally in most conventional gauges, including the 
Newtonian gauge [13]. 



(We imply that tensor perturbations develop entirely 
from primordial metric inhomogencitics, i.e., there are 
no primordial "isocurvature" tensor excitations.) Sub- 
horizon oscillations of the Fourier modes are mapped to 
small-scale discontinuous features of the Green's function 
at the horizon \x\ = r. Only when the phase shift is ab- 
sent [ipo of eq. (35) is zero] the corresponding, step-like, 
discontinuity [eq. (36)] is consistent with condition (4). 

The amplitude suppression factor Aq is determined by 
the Green's function jump at the particle horizon. In 
the order linear in R v = p v /p, when neutrino stress is 
calculated ignoring its feedback to the metric, we find 

A Q = 1-^R V . (5) 

For three neutrino flavors {R v ~ 0.405 for N vc s ~ 
3.04 [14-16]) the agreement of eq. (5) with a numerical 
value 4 CMBFAST ) ~ 0.76 is 2%. 

Fourier transformation of the Green's function, which 
we derive in the 0(R V ) order, gives an analytic expression 
for the evolution of tensor modes in the radiation era. 
Good match of corresponding equations (34,45) with the 
full numerical CMBFAST [17] calculations is demonstrated 
on the right panel of Fig. 1. 

The amplitude and phase of gravity wave oscillations 
are affected by changes in the Hubble expansion dur- 
ing the subsequent cosmological epochs. On the smallest 
scales (k/Ti — > oo) the amplitude and phase shift evolve 
adiabatically (infinitesimally over an oscillation period). 
Then the amplitude decays as 1/a while the phase shift 
remains frozen. Hence, 

h{r, k -► oo) = A{t) sin if /if = A(t) /i (rad ' 0) (tp), (6) 

where A(t) is given by eq. (52). In real space, this 
small-scale behavior corresponds to a step discontinuity 
of Green's function at the particle horizon, 

/i (disc) (r,a;) = 9(t - \x\) = A{t) /i< rad '°> (r, x). (7) 
2r 

The modes entering the horizon in the matter era (k <C 
fe eq ) encounter negligible anisotropic stress. For them 
during matter domination 

h (mat,0) = Sj^/tp = 3 (sm^/V 3 - cos tp/ f 2 ) . (8) 
The respective Green's function is 

ft (mat ' 0) = £ (l - ^) 9{r - \x\). (9) 

Real space results (7) and (9) suggest a compact an- 
alytic expression that correctly describes the tensor evo- 
lution in both the small and large scale limits and is 
qualitatively acceptable in the intermediate range: 

h ~ A(t) h^ Iad ' 0) + [1 - A(t)} M mat >°). (10) 

At any r, the first term guarantees the proper k — > oo 
behavior (6), while on superhorizon scales h(k — > 0) = 1. 
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(The real space equivalent of the latter condition is 
J h(x) dx = 1.) In the matter era, in agreement with 
further accurate analysis, the phase shift of subhorizon 
modes varies with k but is constant in time, eq. (57). Ex- 
pression (10) reduces to ft,( mat '°) for a mode that enters 
the horizon in the matter era, when A(t) decays. This 
simple parameterization matches well the numerical ten- 
sor evolution to the CMB last scattering and beyond, 
Figs. 2 and 3. Yet, despite the close match, formula (10) 
is only qualitative for a finite k. A systematic treatment, 
presented further, reveals important effects eluding this 
minimal description. 

As an alternative to eq. (10), a very popular parame- 
terization of tensor evolution by Turner, White and Lid- 
sey [18] considers an ansatz h TWL = T(k/k cq ) h( mat >°\ 
Here, a "transfer function" is T(y) = [1 + I. My + 
2.5?/ 2 ] 1 / 2 , with its last two coefficients fixed by a numer- 
ical fit at z = 0. This parameterization is rather mis- 
leading. Most importantly, the ansatz T(fc//c oq ) /i( mat '°) 
is motivated by an erroneous statement that "for r 3> 
T eq , the temporal behavior of all modes is given by 
3ji(fcr)/fcT". On the contrary, both subhorizon solu- 
tions of the wave equation in the matter era, ft,( mat <°) — > 
—3cos(p/(p 2 and the complimentary n — > 3sintp/tp 2 , 
eq. (80), decay as l/a. A mode with k > k eq for r 3> T cq 
is described by their linear combination, which is deter- 
mined by the evolution before matter domination. In 
fact, only the complimentary solution contributes when 
k 2> k eq , eq. (6). This ironic feature and the importance 
of the phase shift for the CMB signatures of gravity waves 
was emphasized by Ng and Spcliotopoulos [19]. 

Other shortcomings of the TWL transfer function 
are disregard of the neutrino-induced suppression of 
the modes with k > fc oq and of the dynamical effects 
caused by the residual radiation density around the CMB 
last scattering (p la d/p | 7 doc ~ 0.25 for the WMAP pa- 
rameters [20]). The second problem was alleviated by 
Wang's [21] time-dependent transfer functions which con- 
tain two additional fitting parameters. Nevertheless, his 
parameterizations continued to evolve to the incorrect 
cos ip/<p 2 form even for the shortest wavelengths. 

We obtain the functional form of the tensor evolu- 
tion at a finite k by analytically solving the dynamical 
equations under controlled approximations. In addition 
to leading naturally to adequate simple analytical de- 
scription, derivation, as opposed to fitting, also offers 
better means of identifying the physics responsible for 
evolution features. For earlier works that applied differ- 
ent methods but were guided by the same intentions see 
Refs. [19, 22, 23]. 

The modes that enter the horizon before equality 
(fcr cq > 1) offer a surprise. Contrary to intuitive expec- 
tations, the phase shift of their subhorizon oscillations 
grows after the horizon entry throughout the rest of the 
radiation epoch. The growth is due to the increasing role 
of nonrelativistic matter in the Hubble expansion. When 
matter becomes dominant, the phase shift saturates at a 
value (74). 



We analyze the subhorizon evolution with the WKB 
expansion in the orders of k/Ti.. To start with an order 
above the adiabatic approximation, we factor out the l/a 
decay by considering a variable v cx ah, eq. (49) [24] . Un- 
der the horizon, the amplitude of v oscillations freezes. 
However, their phase shift grows during radiation dom- 
ination as hiy + const, where r e = T cq /(\/2 — 1) ~ 
260 Mpc is given by eq. (47). The oscillation amplitude 
and the additive constant in the phase are set by the 
mode evolution during the horizon entry, described next. 

Matching the ncxt-to-adiabatic WKB order to the su- 
pcrhorizon primordial fluctuations is rather subtle. For 
initial oscillations the WKB approximation breaks down 
and v WKB diverges. Ng and Speliotopoulos [19] ap- 
proached the problem by changing the evolution vari- 
able to Iiit, now varying from — oo. They matched the 
growing mode in the "forbidden" supcrhorizon region to 
the "allowed" subhorizon one by applying the standard 
WKB matching procedure at a turning point. Their re- 
sult reproduced accurately the phase of subhorizon oscil- 
lations but significantly underestimated the amplitude. 
Pritchard and Kamionkowski [23] stretched the WKB 
approach by retaining the decaying mode in the forbid- 
den region. The obtained in terms of special functions 
amplitude remained in poor agreement with numerical 
calculations. 

We consider the WKB solution for v(t) only on sub- 
horizon scales. To normalize its amplitude and fix the 
phase, during the horizon entry we, instead, treat the 
matter correction to the Hubble expansion perturbativcly 
in /O m at/p- Both the WKB and perturbative approxima- 
tions are valid and both solutions can be matched when 
a mode has entered the horizon but the matter density is 
still negligible. The result, in terms of elementary func- 
tions, 

h = A(t) sin[<p - A<p(t)]/<p (11) 

with the amplitude suppression (72) and phase shift (70) 
is in excellent agreement with CMBFAST computations 
(Fig. 4) . As a comparison, the matching procedure of [23] 
reproduces the amplitude of subhorizon oscillations for 
k = 0.1677 Mpc -1 to 13% [23]. For the same k and same 
cosmological parameters, our formula achieves 1%. 

The evolution of the modes that enter the horizon af- 
ter equality (fcr cq < 1) is likewise affected by the sub- 
dominant radiation density. The effects in this case 
are less dramatic. In the linear in kr eq order the im- 
pact of the residual radiation is fully accounted for by a 
shift r — > t + r e in the radiation-free transfer function: 
h{r) w fr( mat >°)(T + r e ). This result is obtained by ob- 
serving that for r < r cq the large-scale scale modes are 
essentially frozen and that for r r oq the background 
expands as a(r) sw a( mat '°)(r + r e ). The advance of the 
tensor evolution on the large scales by the constant con- 
formal time increment T e is consistent with numerical 
computations (Fig. 5). 
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Tensor metric perturbations were connected to the 
CMB temperature and polarization spectra in the origi- 
nal works [25-28] , a full systematic approach was devel- 
oped in [29], and a review of the physical picture and 
typical analytic approximations given recently in [23]. 
We do not elaborate this connection in the present ar- 
ticle. Yet, we remember that generation of CMB polar- 
ization and its curl- type B component, clean from lin- 
ear scalar contribution [5, 6], requires photon scatter- 
ing on free electrons, to polarize photons, but moder- 
ate optical depth, for the polarized photons to reach the 
Earth. These conditions are met inside the surface of 
CMB last scattering (z ~ 1090 ± 100) and after reion- 
ization (z < 20 — 6). The tensors also contribute to 
CMB temperature anisotropy, primarily, through the in- 
tegrated Sachs- Wolfe mechanism [25], induced by hij af- 
ter the last scattering (z < 1090). 

Among other results of this paper, we note sim- 
ple equation (15) that describes fully general-relativistic 
transport of radiation of decoupled ultrarelativistic par- 
ticles. The radiation intensity is quantified by a vari- 
able 7(aj^,ni), uniquely specifying the particle energy- 
momentum tensor with eq. (13). An advantage of the 
suggested intensity definition /, given in Appendix A, 
is that I transport is affected only by spatial gradients 
of the metric and not by metric temporal change. As 
a result, the intensity and its statistical distribution are 
time-independent on superhorizon scales in any gauge. 
This fully nonlinear conservation does not require adia- 
baticity of cosmological perturbations. If particles have 
ever been in thermal equilibrium and decoupled on super- 
horizon scales, a superhorizon value of their intensity I 
is simply related to metric perturbation on a hypersur- 
face of uniform particle density, eq. (A20). In the locally 
Minkowski metric / reduces to the regular particle inten- 
sity dE/(dVdn fl ). 

The rest of the paper is organized as follows. Sec. II 
reviews the dynamical equations. Sec. Ill solves the cou- 
pled evolution of gravity waves and neutrinos during the 
radiation era. Sec. IV considers the evolution of tensor 
perturbations through the radiation-matter transition to 
the matter era. In Appendix A we describe propagation 
of decoupled particles in an arbitrary metric and in linear 
theory. In Appendix B we confirm by explicit integration 
of the real space dynamics that no amount of free stream- 
ing neutrinos can change the gravity wave phase in the 
radiation era. 

In this paper, space-time coordinates x M = (r, x) cor- 
respond to the perturbed Friedmann-Robertson- Walker 
metric 



compare the analytic formulas with numerical computa- 
tions assuming a flat model with f2 m = 0.3, O^/i 2 = 0.022, 
and h = 0.7. 



II. DYNAMICS 

The linear evolution of gravity waves is described by 
wave equation (1). The gravity waves are driven by neu- 
trino anisotropic stress, which develops after neutrinos 
decouple (z < z„dcc ~ 10 10 ). The driving is dynam- 
ically relevant while neutrinos contribute noticeably to 
the total energy density (z > z cq rts 3200). At those red- 
shifts, all the three standard neutrino generations may be 
assumed ultrarelativistic, given the current cosmological 
bounds on the neutrino masses [20, 30-36]. Then neu- 
trino dynamics is approachable in the full general relativ- 
ity for an arbitrary metric, as discussed in Appendix A 
and summarized by Sec. II A next. 



A. General relativistic free streaming 

We describe the transport of ultrarelativistic neutrinos 
by intensity /, defined in Appendix A as their energy- 
integrated phase space distribution. /(a; M ,7ij) is a func- 
tion of spacetime coordinates and direction of neutrino 
propagation. The direction is specified by the covari- 
ant spatial components rij of a null vector n^ 1 oc dx^ jdr 
which is normalized as Y^=i n f = 1- 

The impact of neutrinos on the metric is determined 
by their energy-momentum tensor Tjf . For a given inten- 
sity I, 




(13) 



eq. (A9), where no(rij) can be found from the null con- 
dition n tll n ll g iiv = 0. Provided the considered perturba- 
tions are superhorizon during neutrino decoupling, the 
initial conditions arc set on a hypcrsurfacc of uniform 
neutrino density by (A20) 



I {u \x, ni ) = 



const 



(14) 



Afterward, the neutrino intensity evolves according to 
transport equation (All) 



, di_ = l 

dx^ 2 



dl_ 

dm 



(15) 



(r) K 



(12) 



The considered intensity I is gauge dependent. 4 How- 
ever, by separating the particle and metric dynamics, the 



with rj^u = diag(— 1, 1, 1, 1). Zero of conformal time r is 
set by the condition r(a — > 0) = in the radiation era. 
The scale factor a = 1 at present. Background curvature, 
if any, is ignored for the considered redshifts. Figures 



4 Similarly to the majority of the traditional cosmological vari- 
ables, e.g. [7], I^x^jUi) defines a "gauge-independent" variable 
in any fixed gauge. 
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presented description achieves, first, compact transport 
equation in the full theory. Second, time independence 
of the intensity for supcrhorizon inhomogencitics, with 
negligible g^^. (The metric, which is partly determined 
by a gauge condition, may nevertheless evolve even on 
superhorizon scales.) In a locally Minkowski metric, / 
reduces to the conventional intensity of neutrino radia- 
tion. 

In linearized theory, the intensity perturbation X = 
SI 1 1 obeys transport equation (A13) 

X + rijJj = 2n l n^n u h liu , i , (16) 

with = (l,Tij) and initial conditions (A21), 



(17) 



The perturbed components of the neutrino energy- 
momentum tensor are determined by X from eqs. (A14)- 
(A17). 



B. Tensor modes 

Tensor linear perturbations over the flat Friedmann- 
Robertson- Walker metric satisfy the conditions hoo = 
hoi = and 



dih tj = 0, 



0. 



(18) 



The tensor perturbations evolve according to a wave 
equation, e.g. [7-9], 



V 2 /i„ = 167rGa 2 S 



'j ■ 



(19) 



where they are sourced by the tensor (transverse) com- 
ponent of anisotropic stress 



3 V*" 



(20) 



Neutrino contribution to anisotropic stress follows 
from eq. (A17) by substituting the explicit linear solution 
for neutrino intensity (A23). When r- ln in this solution is 
set to zero and the above tensor conditions are applied, 
eq. (A23) gives [10] 



X(t,x) — nirijhij(T,x) — / dr' mnjhijfr', x'), (21) 



1 



with x' = x + n(r' — r). 

The impact of photon anisotropic stress on tensor met- 
ric perturbations is not entirely negligible. While local 
photon anisotropy is washed out by Thomson scatter- 
ing before the recombination, the anisotropy develops af- 
ter photons decouple at z 7 dcc ~ 1090. At this redshift, 
radiation constitutes as much as a quarter of the total 
energy density. The generated photon stress affects the 
modes that enter the horizon shortly after the decou- 
pling. Numerical calculations show up to 7% consequent 



suppression of the metric modes with k ~ 0.005 Mpc -1 . 
However, since the photon stress was insignificant during 
decoupling, CMB polarization is almost unsuppressed. 
CMBFAST predicts less than 0.5% suppression of the po- 
larization correlation Cf B . The photon stress also affects 
little CMB temperature - the corresponding suppression 
of the tensor component of Cj T is less than 1%. We 
therefore neglect photon anisotropic stress in the follow- 
ing analysis. Even more so, we can disregard the stress 
of nonrelativistic species (baryons and CDM) because of 
their tiny velocity dispersion. 

Now we consider a single plane wave, not necessarily 
harmonic, and choose the coordinates y and z in a plane 
of constant hij: hij = hij(r,x). Then for tensor pertur- 
bations, by the first of conditions (18), h x % = hi X = and 
Y, x i = Ynx = 0. For the remaining components of 
with i,j € {y,z) we can take the average in eq. (A17) 
over the orientation of (n y ,n z ) under a fixed n x = /i. 
The result is 



dr' / d/x(l-/i 2 ) 2 ^(T>')> (22) 



where x' = x + /i(r' — r). Fourier transformation of this 
expression reproduces eq. (16) of Ref. [10]. 

From now on, we imply a decomposition of a gravity 
wave hij over some polarization basis, e.g., into -I- and x 
components, 



X=l,2 



(23) 



Since the evolution imposed by eqs. (19) and (22) is iden- 
tical for all the polarization components, we consider an 
arbitrary component h\ and drop the subscript A later. 



C. Plane tensor Green's functions 

The evolution of any linear perturbation can be de- 
scribed by superposing plane Green's functions [37, 38] 
h(r,x), which satisfy the initial condition 



h(r ->0,x) = 6(x). 



(24) 



For example, a harmonic plane wave mode h(r, k) e lkx 
with the initial condition h(r — ► 0, k) = 1 is described by 
the Fourier integral 



h(r, k) = 



dx ( 



h(r, x). 



(25) 



Since hAh decays, or "squeezes" [39], on superhorizon 
scales, the initial h need not be specified precisely. It is 
sufficient to require a finite r — > limit of dxh(x), 
to ensure that we work with a non-decaying mode. 

The initial condition (24) has several immediate con- 
sequences. First, by the causality of gravity wave propa- 
gation 



h(r, x) = for |a;| > r. 



(26) 
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Second, by parity conservation during the considered 
evolution, the initially even perturbation (24) remains 
even: h{r,—x) = h(r,x). Third, the k — > limit of 
eq. (25), time-independence of h(r,k — ► 0), and initial 
condition (24) add up to a sum rule 

/oo 
dx h(T, X ) = h{T,k -> 0) = 1 (27) 
-oo 

for all t. 

III. RADIATION ERA 

The magnitude and phase of the tensor modes with 
k 3> fc C q is set by their dynamics in the radiation epoch. 
At that time neutrino anisotropic stress plays an active 
role. The coupled dynamics of gravity waves and neu- 
trinos during radiation domination can be solved analyt- 
ically in real space. Analytic expressions for the radia- 
tion era evolution of Fourier tensor modes then follow by 
Fourier transformation. 



where 9 is the Hcavisidc step function and 

g( ,, x)5 (l^gx%) . (32) 

P>~X 

Equation (31)-(32) is singular at x = arL d X = A*- 
The real-space singularities are treated consistently with 
the Fourier mode approach by the standard formalism 
of generalized functions (see discussion in Sec. IV. C of 
Rcf. [12]). In the following analysis of the tensor sec- 
tor it is sufficient to remember that, first, the d/j, inte- 
gral in eq. (31) should be taken in the Cauchy principal 
value. Second, an equation (\ — a) A = B for any gen- 
eralized functions A(x) and B(\) should be resolved as 
A = B/(x — a) + const S(\ — a), where the last term is the 
Dirac delta function. The delta function prcfactor can be 
determined from the properties imposed on A by the ini- 
tial conditions; for example, from the sum rule (27). 



B. Phase shift 



A. Self-similarity 

The coupled evolution of gravity waves and neutri- 
nos with the specified Green's function initial conditions 
involves no external dynamical scales during radiation 
domination, when TL = 1/r. Then the plane tensor 
Green's function, whose dimension is t _1 , takes a self- 
similar form 



h(r, x) 



Hx) 



x 

X = -• 



(28) 



r t 
The initial condition (24) requires the normalization 

d X h(x) = l, (29) 

as evident from the sum rule (27) . Note that integration 
over x is restricted to x £ [ — 1, 1] because by causal- 
ity (26) h(x) vanishes identically beyond this interval. 

We determine h{x) by solving eqs. (19) and (22). We 
use that h = -{ X h)'/T 2 , h = (x 2 ^)"M V 2 h = h"/ T 3 , 
where primes denote x derivatives, and in the radiation 
era 8nGa 2 = 3/(pr 2 ) by the Friedmann equation. With 
these substitutions, eqs. (19) and (22) give: 



[(x 2 -i)h']' = 

2>R U f T rdTi 



(30) 



d/x(l-^) 



2^2 d [Xl^Xl)} 



lo <i J-i d Xi 

where xi — (x ~ A*)^" + A* anc i — Pu/P- Carefully in- 
tegrating over dn by parts and remembering the causal- 
ity (26) we obtain: 

m 

~2 



0(1-1x1) / d^[K(jM, X )-K(x,n)], 



Without neutrino anisotropic stress (i?„ = 0) , the even 
normalized solution of eq. (31) is 



^ (0) = ^(i-|xD- 



(33) 



Fourier transformation of the corresponding Green's 
function (28) gives the standard neutrinoless radiation 
era result 



sin(A:T) 
kr 



(34) 



Although neutrinos in the radiation era are, in fact, 
among the dominant species, the gravitational impact of 
their anisotropic stress becomes negligible on subhorizon 
scales. The corresponding asymptotic (kr ^> 1) homoge- 
neous solution of mode evolution eq. (1) is 



h(kr) 



A sin(/cr + ip Q ) 
kr 



o ({k T y 



(35) 



We prove that, regardless of an R u value, ipo = 0. 

The amplitude and phase of the transfer function os- 
cillations in kr are fixed fully by the Green's function 
discontinuity at |%| = 1. The fcr — > oo behavior of the 
Fourier modes (35) maps to the following discontinuous 
terms in real space [40, 41]: 



^ (disc) (x) = ^[cos <PO0(1-\X\) 



sin ipo 



ln|l- X 2 |].(36) 



Regardless of the continuous part of h(x), h = for 
|x| > 1, as required by causality, only if sini^o = 0, i.e. 
if the logarithmic term in eq. (36) vanishes identically. 
Appendix B confirms that the explicit solution of eq. (31) 
has this property for any R v . 
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FIG. 1: Evolution of tensor metric perturbations in the radiation era as described by a real-space Green's function, left, and 
Fourier modes, right. For both presentations, the displayed solutions are the neutrinoless, unsourced by any anisotropic stress 
(dashed) and the linear in p v /p for three neutrino flavors (solid). The Fourier modes on the right panel are additionally 
compared with the full Boltzmann CMBFAST calculation (bold dots). 



Since sin 1^0 = 0, either <po = or tpo = 7r. If both 
Aq and <po in eq. (35) change continuously 5 in R v and 
(po = at R v = then ipo = for any R v . 



The delta-function prefactors are determined by observ- 
ing that the spikes control the constant term in hS 1 ' = 
J hy>'dx for x G [ — 1> !]■ This term is to be fixed by h 
normalization (29), by which 



C. First-order solution 

We look for an explicit Green's function h(x) as a series 
in the powers of R v . In the zeroth order h is described by 
a top-hat function (33). The 0(R V ) order is addressed 
next; it will match the full numerical Boltzmann evolu- 
tion to a few percent. Since h = for \x\ > 1, only the 
interval x G [ — 1, 1] is considered. 

Substitution of the zeroth-order h into the right hand 
side of eq. (31) gives for the 0{R V ) correction 



( X 2 -!)/>«' 



(37) 



3i?„ 
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2x 4 -yX 2 + l + (x 2 -l) 2 xln 



i-x 
i + x 



After one integration and division of both sides of the 
obtained equation by % 2 — 1, 



h™' = R v 



(38) 



where % e (—1,1). 

As discussed at the end of Sec. Ill A, the division by 
X 2 — 1 creates delta-function spikes in hS 1 ' 1 at x = ±1- 



(39) 



Integrating eq. (38) and fixing the constant term as 
described we find that 



R v 



1 



h (1) - ^ j6 X 4 -23 X 2 + - + 



(40) 



+ (3 X 4 -10 X 2 + 15)xln^-81n(l-x 2 )| 

for x € [— 1) 1] and = otherwise. 

As anticipated, discontinuity at |x| = 1 is of the 
form (36), where <po = and 



AP= lim 2h^ = -lR v . 



IxHi-o 
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(41) 



The obtained Aq determines the subhorizon radiation era 
solution for tensor modes, eq. (35), as 

h(kr) =(l- 5 -R v + (j£)) ^ + 0((fcr)- 2 ). (42) 

We also note that the leading behavior of the <9((/cr)~ 2 ) 
correction follows from the general result (B7), 



5 If we demand that the amplitude and phase change continuously 
in R v , we should, in principle, allow the amplitude to become 
negative. The following calculations show that in a realistic range 
of R v the amplitude varies moderately and remains positive. 



. sin fcr _ cos fcr rH , . . 
h{kr) = A ik-ro-2- [1 + (!)] . (43) 



kr 



k 2 T 



valid in the radiation era at all orders of the R u expan- 
sion. 
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It is possible to find the entire radiation era evolution of 
a tensor Fourier mode h(kr) in 0{R V ) order. Introducing 
ip = kr, we calculate 



with hy~'(x) given by eq. (40). The result is 

h W (<p) = ??-( A« sin ip - cos 
V v 

where 



(44) 



(45) 



=-l-j* + ^ + ^ CS(2v») - 4 - 81(20, 
* (1) = i + ^ + £ffl(2*0+4(£-£) Ci(2^. 



Here, Cix = J* — - dy = — hix — 7 + cix and 



Six = Jq ^Y*-dy = § + six. Explicit verification with 
Mathcmatica confirms that in the 0(R V ) order the ob- 
tained /i^ - 1 + M 1 ' satisfies Weinberg's integro-differcntial 
equation for a tensor Fourier mode, eq. (21) in Ref. [10]. 

Figure 1 shows the neutrinoless (dashed) and the linear 
in R v (solid) gravity wave solutions in real and Fourier 
presentations. Here, R v is 0.405, corresponding to the 
three standard neutrinos. The Fourier modes are addi- 
tionally compared with a full Boltzmann calculation with 
CMBFAST [17] (dots). As seen, the 0(R V ) analytical so- 
lution provides a reasonable accuracy for all kr. 

For the standard R v value, the 0(R V ) tensor amplitude 
suppression factor is A = 1 — ~ 0.77. It is in ex- 
cellent agreement with the CMBFAST calculation, giving 
^(Cmbfast) ^ 0.76 during radiation domination (red- 
shift z = 10 7 , the same result at z = 10 8 ), and with nu- 
merical solutions of Weinberg's intcgro-diffcrential equa- 
tion in Fourier space, _ i 4( Weinb ) ~ 0.8O [10] and Aq PK - ) ~ 
0.81 [23]. 



IV. TO THE MATTER ERA 

Most of the detectable gravity waves effects are gen- 
erated after the radiation era. The modes probed by 
interferometers and pulsar timing enter the horizon in 
the early radiation era, and their subsequent evolution 
is described excellently by the adiabatic approximation 
h tx sin(fcr)/a. The modes probed by CMB anisotropy 
and polarization enter when radiation and matter densi- 
ties are comparable, and the post-radiation era evolution 
of these modes is more subtle. Note that their CMB im- 
prints are left primarily at redshifts at which dark energy 
plays little role in the background expansion. 



A. Minimal analytic solution 

For the studied redshifts, when dark energy density can 
be assumed small, we describe the Hubble expansion by 



a radiation-matter model with p = [p 7) o/(l — Ru)]/a A + 
Pm,o/a< 3 - The corresponding scale factor a, governed by 
the Friedmann equation, evolves as 



2f 



T 2 ) £ 



eq ■ 



(46) 



Here, a cq = [p 7 .o/(l — Rv)}/ Pm,o is the value of a at 
radiation-matter equality, and f is conformal time in 
units of characteristic conformal time of equality 



t/t c 



(47) 



(>/2-l)r e 



Note that by eq. (46) r(a oq ) = t, 

First, we address the modes that enter the horizon long 
before equality (k 3> t~^) although after neutrino decou- 
pling. After a mode has become subhorizon (k 3> Tl), we 
ignore the anisotropic stress. Then the mode evolves ac- 
cording to a homogeneous wave equation 



h + 2Hh + k 2 h = 0. 



(48) 



The Hubble rcdshift 2Hh of the graviton energy is ac- 



counted for by substitution h = 



k 2 - - \v = 0. 

a , 



diere now [24] 



(49) 



In the small-scale limit k/Tt — > 00, the term a/a ~ H 2 
can be dropped. In this (adiabatic) approximation v os- 
cillates harmonically in tp = kr with constant amplitude. 
Correspondingly, the amplitude of h oscillations decays 
as 1/a. 

We normalize the amplitude and fix the phase shift of 
the considered short-scale modes by subhorizon radiation 
era solution (42). As a result, 



h(r,k 



A(r)hW>(<p), 



where 



sirup 



(50) 



(51) 



is the neutrinoless radiation era solution, and A(r) oc tp/a 
equals 



A(r) 



1 



R v 



(52) 



The last formula neglects the 0(R^) correction in eq. (42) 
and uses a(r) form (46). We note that the Fourier short- 
scale evolution (50) maps onto the following Green's func- 
tion discontinuity, c.f. eqs. (35)-(36): 

/i( disc )(r,x) = ^-6(T-\x\)=A(r)h ( - ta - d ' 0) (T,x). (53) 
2r 

Next, we consider the opposite, large-scale limit (k <C 
tJ^ 1 ), in which mode evolution starts only in the mat- 
ter era. Then anisotropic stress of relativistic species is 
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negligible and, similarly to the radiation era, the tensor 
dynamics is scale-free. Solving eq. (19), where now X = 
and H = 2/r, with a normalized Green's function of the 
form (28), we obtain 

^(mat.O) = 3 L _ X \ fl , (M) 

4r \ t z J 

Fourier transforming, we arrive at a well-known tensor 
transfer function in a matter dominated universe 

^ = 3^)^ = 3 (5^-^). (55) 

Finally, we combine the Green's function discontinu- 
ity (53), describing the short modes, with the matter era 
evolution (54) on the larger scales by a "minimal" ana- 
lytic formula. Namely, we look for an expression of the 
form (a + bx 2 ) 9{r — that is normalized (J hdx = 1) 
and jumps to A(t)/(2t) at the horizon = r. Such an 
expression is unique and is given by 

h ~ A(t) h^ dfi) + [1 - A(t)} ^ mat >°) . (56) 

If we view formula (56) as a transfer function h(r, k) 
in Fourier space, with /i( rad <°) (51) and /i( mat <°) (55), we 
observe that h approaches unity for cither r — > or k — > 
0. This formula also reproduces the correct, in 0(R U ), 
asymptotic behavior for k — > oo at any fixed r. 

Figure 2 compares the approximation (56) (solid curve) 
with numerical CMBFAST [17] calculation (bold dots) at 
CMB last scattering, z = 1090. The heights of all the 
oscillation peaks match within 11% and the phase within 
3% of 7r. Neither the radiation-era transfer function 
A(t) /i( rad > ' , corrected for the neutrino and adiabatic 
damping (dashed) nor the popular matter-era transfer 
function M mat, °' (dash-dotted) alone can describe the 
tensor evolution to CMB last scattering satisfactorily for 
all k. 

The good agreement between minimal analytic for- 
mula (56) and numerical calculations is preserved at later 
times, which are relevant for tensor signatures in ISW dis- 
tortion of CMB temperature and large-scale polarization 
from reionization (reionization bump) . A comparison for 
a redshift z = 109 is given in Fig. 3. 

We highlight that approximation (56) preserves the 
temporal phase shift of mode oscillations after radiation- 
matter equality. On subhorizon scales, this formula re- 
duces for t ^> T eq to 

h (sub, mat) 2fcT e(! ~ | jk) Sin <g - 3 COS tp 

If 2 

As the original formula (56), asymptotic expression (57) 
has the correct leading-order behavior in both kr e — > 
and kr e — > oo limits. However, the actual subleading 
corrections to these limits differ from the terms of eq. (57) 
in magnitude and, when kr e — > oo, even in k scaling. The 
evolution of a tensor mode with a finite k is studied next 
quantitatively 
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FIG. 2: Tensor modes at the redshift of CMB last scattering, 
z = 1090. The plots are generated with: "minimal" analyti- 
cal formula (56) (solid) , numerical CMBFAST calculation (bold 
dots), the neutrino- and adiabatically-damped radiation era 
solution (dashed), and the matter era solution (dash-dotted). 
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FIG. 3: The same calculations as in Fig. 2 but for a ten times 
smaller redshift. "Minimal" analytical expression (56) (solid) 
continues to provide good accuracy. 



B. Entering before equality (fcr cq > 1) 

First, to determine the role of the radiation-matter 
transition in mode evolution, we ignore neutrino 
anisotropic stress. The derived formula will ultimately 
be corrected for it. Without anisotropic stress, an exact 
formal expression for tensor evolution in the radiation- 
matter background can be given in terms of radial prolate 
spheroidal functions [22, 42]. Yet, a cumbersome form of 
this expression (an infinite series of Hankcl functions with 
coefficients involving more infinite series), calls for more 
insightful and practical description. 

We continue to parameterize tensor mode evolution 
by the dimcnsionlcss oscillation phase tp = kr and use 
the perturbation variable v oc ah. We normalize v by 
v(tp) = (p + o(tp) as tp — > 0. Then the initial condition 
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h(r — > 0) = 1 and explicit a(r) solution (46) fix ah/v to 

h= ^ v= !L_ . (58) 

atp e (1 + if) if 

In the first of these equalities and later 

tp e = kr e . (59) 

Source-free dynamical equation (49) for v(tp) in the 
radiation-matter background reads 



+ 1- 



1 



We 



y 2 



v = 0. 



(60) 



Here and below primes denote derivatives with respect 
to tp. In a purely radiation background, in which a = 
and v" + v = 0, the solution is ?/ rad > ) = simp. To 
account for massive matter deep in the radiation era 
(ip -c tp e ) we treat (d/a)v in eq. (49) perturbatively 
approximating it by (d/ a)v^ Tad '°\ Then, neglecting tp 2 
against tptp e , we have 



v + v 



sin 93 

We ' 



(61) 



A substitution v = A(tp) sin tp leads to a first-order dif- 
ferential equation for A' (ip). That equation has a unique 
regular at tp = solution 

A' = - 1 2 pn(2v>) + 7 - ci(2 V )] , (62) 
2iy9 e sin 93 

where 7 = 0.577216... is the Euler constant and 
ci.T = — J cosydy/y is the cosine integral. Integrat- 
ing eq. (62), we obtain: 



w (rad) _ A(tp) sin tp - B(tp) cos tp, 



where 



.4 
B 



ltp e 



si{2tp) 



— \H2ip) + 7 - ci(2<^)] 
2ip e 



(63) 

(64) 
(65) 



and six = — sin ydy/y is the sine integral. 

Well under the horizon (tp 3> 1), the sine and cosine 
integrals in expressions (64)- (65) tend to zero. Then 



/rad,\ 
y\ sub / ; 



7T \ 1 

- — sin tp - — [ln(2tp) + 7] cos tp. (66) 
Atp e J 2tp e 



We see that at a large but finite tp e = kr e the oscillation 
phase is shifted from the sin tp oscillations of the adia- 
batic limit. The shift is caused by nonrelativistic matter 
contributing to the background expansion. As the mat- 
ter role increases in time, the shift grows logarithmically 
even after the tensor mode has become subhorizon. 



We can easily extend the mode subhorizon evolution 
beyond the radiation era using the WKB approach 6 . To 
that end, we combine the prcfactors A and B of eq. (63) 
into a varying complex amplitude C = A — iB: 



Im(C. 



(67) 



C(tp) satisfies a second-order differential equation, follow- 
ing from eq. (60). After the horizon entry C(tp) changes 
slowly, hence, we neglect C"(tp) in that equation. The 
remaining terms are 



c 



c. 



(68) 



2tptp e + tp 2 

Integrating (68) and setting the integration constant to 
match the mode (67) to radiation era solution (66), we 
find 



w (sub) = 1 



with 



^ 2^ e 



4ip e 



-In 



sin(( / 9 — Atp), 



1 

itp e 



2<p 



7T- +7 



(69) 



(70) 



Finally, we incorporate the damping by neutrinos dur- 
ing the horizon entry in the radiation era. We treat the 
correction caused by neutrino anysotropic stress as per- 
turbation which is independent of the perturbation in- 
duced by the matter term in eq. (61). In this approx- 
imation, neglecting interference of the two corrections, 
neutrinos additionally suppress the oscillation amplitude 
by the factor in eq. (42) but they do not affect the oscil- 
lation phase. Adjusting the ?/ sub ) amplitude in eq. (69) 
and substituting the result into eq. (58), we thus have: 

, A(r,k) 



sm[tp - Atp(r, k)], 



where on subhorizon scales 

^(sub) ^ 



(71) 



(72) 



and Atp is given by eq. (70). 

Deep in the matter era > tp e ), the obtained A( sub '' 
decays inversely with time as 



^(sub, mat) 



2(^(1 -§i?„) + § 
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while the A<£>( sub ' growth saturates at 

,(sub,mat) _ ln(4v? e ) +7 



Atp^ 



2tp e 



(73) 



(74) 



Since tp e = kr e , the phase shift of tensor oscillations at 
fixed time and k — > 00 approaches its asymptotic zero 



6 Since the adiabatic decay /i ~ 1/a on subhorizon scales is already 
accounted for by the v prefactor, e.g. eq. (58), the presented first- 
order WKB calculation of v is equivalent to the second-order 
WKB analysis of the tensor metric perturbation h. 
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FIG. 4: The sixth oscillation of a tensor mode with k = 
5t~^ as: next-to-adiabatic analytical solution (71)-(72)- 
(70) (solid), CMBFAST output (dots), and the neutrino- and 
adiabatically-damped radiation era solution (dashed). The 
top figure is for the p v — model. Both the matter-induced 
amplitude enhancement and phase shift are seen well repro- 
duced analytically. The bottom figure describes the stan- 
dard 3-neutrino scenario. The appearing small discrepancy is 
caused by the interference of the neutrino and matter correc- 
tions, as discussed in the main text. 



value as In k/k, more slowly than could be expected from 
naive Taylor expansion. We recall the origin of this scal- 
ing in the impact of matter density, which fraction in- 
creases in time and induces logarithmic growth of the 
tensor phase shift before radiation-matter equality. 

The accuracy of our result (71)-(72)-(70) against nu- 
merical calculations is seen from Figs. 4, which capture 
the sixth oscillation cycle of a tensor mode with fcr cq = 5. 
The plotted are: the neutrino- and adiabatically-damped 
radiation era solution (dashed), our next order analyti- 
cal prediction (solid), and the full CMBFAST calculation 
(dots). The top Figure describes the p„ = case and 
shows the excellent accuracy of the next-to-adiabatic an- 
alytical predictions for both the amplitude enhancement 
and phase shift induced by matter. The bottom Fig- 
ure refers to the standard 3-neutrino scenario. Then, 
in addition to the neutrino and matter effects, their un- 
accounted interference slightly increases the amplitude 
and the phase shift. The higher amplitude is expected 
from dilution of neutrinos by matter. The larger phase 
shift may be understood as follows: in the radiation 



era according to eq. (43) neutrino anisotropic stress in- 
duces a decaying shift of oscillations toward higher tp. 
When the neutrino source is cut off by matter domina- 
tion at a finite tp ~ tp e , the oscillations remain shifted by 
&„<p = 0{R v /<p e ). 



C. Entering after equality (kr cq < 1) 

Now we address the modes that start to evolve after ra- 
diation becomes subdominant to nonrclativistic matter. 
We determine the tensor evolution by retaining only the 
leading, linear in <p e = fcr e , corrections due to the resid- 
ual effects of radiation. We could treat the radiation- 
induced corrections perturbatively, similarly to treating 
the subdominant matter effects in the previous subsec- 
tion. However, the derivation becomes more elegant and 
the result is easier to interpret if, instead, we account for 
the radiation energy by a simple change of variables, as 
described next. 

We start by rewriting the scale factor evolution (46) as 



a c 



[(T + Te) 5 



/rl 



In terms of the evolution variable ip = fcr, 



a = 



where 



(p = fc(r + T e ) = tp + ip e 



(75) 



(76) 



(77) 



Next, in eq. (76) we drop the last tp\ term and in the 
gravitational wave eq. (1) ignore the anisotropic stress 
of radiation, both being O(tpl) corrections. Now, with 
a(r) (x p 2 , we observe that all the remaining contribution 
of radiation to eq. (1) can be accounted for by considering 
a matter-only scenario but replacing tp by the shifted 
evolution variable tp, whose range is (tp e ,+oo). In other 
words, a solution of full eq. (1) can be written as 



h(r,k) = U mat ' \<p) + O(cpl), 



(78) 



where ft,( mat '°) is a matter-only solution with so far un- 
specified initial conditions. 

Although the omission of the last tp 2 in eq. (76) has 
a big impact on H. when t < r e , the considered large- 
wavelength perturbation is then frozen, whether we apply 
the true or modified Hubble expansion. Therefore, as can 
be easily verified, the p\ omission leads to a minor 0{tp 2 e ) 
change in the mode evolution. 

In view of the last remark, we continue to impose the 
initial conditions h = 1 and h' = at ip = 0, i.e., at 
tp = tp e . These conditions fix the constants C\ and C2 in 
the general solution 



h(r, fc)«C , 1 /i( mat ' )(^) + C 2 n(<^), 



(79) 



where /i( mat .°) is the regular matter era solution (55) and 
n is a complimentary linearly independent solution. To 
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— Shifted h (mat ' 0) 




FIG. 5: Evolution of a mode with k = 0.2 r,," 1 as: the shifted 
r — > T + T e matter era solution (solid), the regular matter era 
solution (dash-dotted), and CMBFAST integration (dots). 



APPENDIX A: FREE STREAMING IN 
GENERAL RELATIVITY 

1. Full theory 

Classical particles can be described by a phase space 
distribution / over spacctime coordinates (r, x l ) and par- 
ticle canonical momenta Pi. We consider Pi, as opposed 
to more conventional [43-45] comoving proper momenta, 
as the primary phase space coordinates. At a minor ex- 
pense of two extra terms in the linearized equations for 
the energy-momentum tensor, eqs. (A14) and (A17), our 
choice leads, first, to a simpler general-relativistic equa- 
tion of / evolution. Second, to time independence of / 
for free particles on superhorizon scales [46]. 

The phase space distribution of free streaming particles 
evolves according to the collisionlcss Boltzmann equation 



be specific, we take 
n((p) = 



cos ip sin (p 



(80) 



Then /i( mat >°) and n are comparable on subhorizon scales. 
Given that /i( mat < )(^ e ) = 1 + 0(<p 2 e ) and n(ip e ) = 
0(tp~ 3 ), the initial conditions fix the weights in eq. (79) 
as 

d = l + 0(<£), C 2 = 0(v 5 e ). (81) 
We conclude that up to 0((p 2 ) corrections 



h(T, k) PS /l( mat <°)(<^) 



sin ip cos ip 



(p '■ 



(P 1 



(82) 



This result can also be derived by a regular pcrturbative 
in ip e analysis. 

We thus found that in 0(ip e ) order the residual radia- 
tion density pushes the tensor modes ahead of the zeroth- 
order solution /i( mat, °)(fcr) by a constant conformal time 
increment r e ~ 260 Mpc. Fig. 5 shows good agreement 
between eq. (82) (solid) and CMBFAST integration (bold 
dots) for a mode with k = 0.2 r" 1 . 

In the subhorizon limit, eq. (82) reduces to 



^(mat.sub) _ 3 cog ^ + lpe )/^ 



(83) 



The amplitude suppression factor A and the phase shift 
Aip of our earlier parameterization h = A(t) sin[<£> — 
A<£>(t)]/V are now equal 

^(mat^ub) PS -, A(^ mat ' Sub ) PS — — ip e , (84) 

ip 2 

These values for the large wavelengths can be compared 
with results (73) and (74), derived for the scales that en- 
ter the horizon before equality. If desired, these two sets 
of asymptotic formulas can be joined by various fitting 
ansatzs for a more accurate interpolation than provided 
by our minimal formula (57). 



f,dx^dl,dP L dl = 

1 Sfi At BP, 



(Al) 



The canonical momentum of a particle that has mass m 
and moves along a worldline x M (r) = (t,x 1, (t)) is 



Pi = mg^dx^/i-ds 2 ) 1 ' 2 , 



(A2) 



where ds 2 = g fiV dx^dx v . Natural definitions Po = 
mgo^dx^/i-ds 2 ) 1 / 2 , P^ = g^ v P v and the calculation 
of dPi/dr from the geodesic equation give 



dx l 



P l 



dr 



2P° 



(A3) 



The value of Po at every phase space point {x % , Pi) follows 
from the identity 



(A4) 



We specify the direction of particle propagation by the 
spatial components of 



„ = where P 2 = £ P 2 



(A5) 



This metric- independent, therefore non-covariant, defini- 
tion of separates the infinitesimal transport of parti- 
cles from the temporal change of metric. Note that by 
eq. (A5) 



(A6) 



We also consider = g^ v n v = P M /P. 

The energy-momentum tensor of the described parti- 
cles equals 



TH = 



d 3p. pil p v 
i/=g~ P° 



(A7) 
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Wc substitute P u 



n^P and P» 



n^P, and define 



particle intensity I in a direction as 

/>oo 

7(s",n<) = / P 3 dPf{x^n t P). 
Jo 



Then 



d 2 ni n^riv 



i. 



(A8) 



(A9) 



The dynamics of the intensity of ultra- relativistic decou- 
pled particles is governed by a simple closed equation, 
derived next. 

Only two of the three components of rii, constrained 
by condition (A6), are independent. Let /j, a = (/Ji, ^2) be 
any two independent variables parameterizing rij. Then 
partial differentiation with respect to rij = Pi/ P is natu- 
rally defined as a differential operator 



d 

drii 



a=l,2 



d 



(A10) 



Integrating the Boltzmann equation (A1)-(A3) over 
P 3 dP, applying the above definition of d/drii, and re- 
membering the condition of ultra-relativity g^n^n^ = 0, 
we obtain 



01 



1 



n h ~- — = -n ll n u g l _ ivl 
dx» 2 yM ' 



81 

dn. 



^ii-^r . (Ail) 



Equation (All) is the main result of this Appendix. 

As this dynamical equation shows, for superhorizon 
perturbations the considered intensity of decoupled par- 
ticles is time independent. Indeed, when the spatial gra- 
dients dl/dx 1 and <7 M „,i are negligible, eq. (All) becomes 
1 = 0. In the opposite extreme, for the locally Minkowski 
metric, I reduces to the conventional particle intensity 
dE/{dVd 2 n). 



1. Linear theory 



For a homogeneous and isotropic distribution of parti- 
cles with energy density — Tq = p(r) we see from eq. (A9) 
that I = pa 4 /(Aw). We describe linear inhomogeneities 
over this background with relative intensity perturba- 
tion 2(x M , ni), 



(A12) 



When the particles are ultra-relativistic and their 
number is conserved (pa is constant), linearization of 
eq. (All) gives: 



2 + n{L i = 2n i h* J 'h v h 



(A13) 



where = (l,?ij). The linearly perturbed 

components (A9) arc straightforwardly calculated from 



rifiUvg^ = and (m) n = 0, (ninj) n = %/3, 
() n stands for J d 2 n i /47r: 



J n — 



P (l + (2)n - 2/l) 



(riil)n, 



Here, h = \ ^ l h l 



Ei = p 



4 - 



where 

(A14) 
(A15) 
(A16) 

(A17) 



is anisotropic stress, and in the last formula hij = hij — 
Sijh is the traceless part of hij. 

If the energy distribution of particles that stream 
in a direction n, is thermal and the metric is locally 
Minkowski then 2 is proportional to the particle temper- 
ature perturbation: 2 = 4AT(rii)/T. Under a change of 
coordinates x' 1 — > x M = x M -I 
X transforms as 



(gauge transformation), 



I(^,n t )-I(^,n t ) 



-4riin M a^. 



(A18) 



Similarly to the multipoles of AT(rii)/T, the mul- 
tipolcs of intensity perturbation X(rij) can be used to 
decompose the transport eq. (A13) into a computer- 
friendly Boltzmann hierarchy of equations for photons 
(with scattering terms and polarization added), neutri- 
nos, and other cosmological species. These equations for 
scalar perturbations in the Newtonian gauge are given 
in Sec. Ill and Appendix A of [46] ([46] uses D = §2). 
The structure of these equations is considerably simpler 
than that of the traditional description with AT(rii)/T 
multipoles: The dynamical matter perturbations are no 
longer driven by time derivatives of gravitational po- 
tentials, which are constrained non-dynamically by the 
same matter perturbations. Moreover, unlike the situ- 
ation with the traditional approach, the change of the 
perturbation value during the horizon entry is controlled 
only by physical interactions and is not caused by gauge 
artifacts [46]. 



3. Initial conditions 

We specify the initial conditions by following Rcf. [10]. 
We suppose that the considered species (later, neutrinos) 
decouple at time t^ cc when the Hubble scale is much 
smaller than the spatial scale of the probed perturba- 
tions. Then shortly after Td oc there exists a hypersurface 
with the following property. Every observer on this hy- 
persurface who moves normally to it would detect the 
same isotropic distribution of the proper neutrino mo- 
mentum p. This is a hypersurface of constant tempera- 
ture when the particle distribution is thermal. Although 
the decoupled relic neutrinos are not exactly thermal, the 
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described conditions should still be met on a hypersur- 
face (u) of uniform neutrino energy density. We choose 
u for an initial constant time hypersurfacc. 

An observer who moves normally to a constant time 
hypersurface measures proper momentum p which is re- 
lated to the canonical momentum Pi as 



The trivial integration of the full time derivative gives 



1 



dr' h^Luir',^). 



(A23) 



p 2 = ; //'•'/',/', (A19) 

(iS)gij j g ^ ne mv erse of the 3-tensor gr«). In terms of the 
variables (A5), p 2 — ( 3 ^ IJ mrij P 2 . Therefore, if the hy- 
persurface of the uniform distribution of the proper mo- 
mentum is taken for a constant time hypersurface then, 
according to eq. (A8), the neutrino intensity varies as 



const 



In linear theory and the metric (12), ^g 1 -* 



(A20) 



hij). Then comparing eq. (A20) and (A12) and remem- 
bering (A6) we conclude that on supcrhorizon scales the 
intensity perturbation equals 



This solution is applicable to any (scalar, vector, or ten- 
sor) type of linear perturbation of ultrarelativistic species 
which decouple on superhorizon scales. 



APPENDIX B: FORMAL PROOF OF lim <p 

k — >oo 



Sec. Ill B of the main text argues that subhorizon 
tensor oscillations in the radiation era 
Aosiii(ip + ipo)/ip, ip = kr are consistent with causal- 
ity if only ipo — 0. The unshifted oscillations match to a 
finite, step discontinuity of the real space tensor Green's 
function at the particle horizon \ — ±1: 



^(dis 



) (x) = ^^(i-W)- 



(Bl) 



2(») 



(A21) 



Linear solution 



Sourced transport equation (A13) with initial condi- 
tion (A21) at r = r- m is solved by 



-I(t, x,n) = ninjhlj\T m> x- m ) + (A22) 



+ / dT l n l h f *h u h IJ , u , l (T',x'), 



Vice versa, a Green's function discontinuity of this type 
implies that ip = 0. 

In this Appendix we verify that the solution of the real 
space equation (31) for h(x) indeed has the form (Bl). 
We do not assume that R v is a small parameter. To be 
specific, we consider the h(x) discontinuity at x = — !■ 

We take h(x) of the general form (36) and first show 
that the right hand side of eq. (31) should remain finite 
as x ~ + —1 + 0. For the second term in the integral on 
the right hand side of eq. (31), 7 



lim 

x^-i+o 



d^K(x,^) = 0. 



(B2) 



where Xi n = x + n(Tj n — r) and x 1 = x + n(r' — r). 
At the lower limit of the integral, should, strictly 
speaking, be evaluated on the uniform density hypersur- 
face on which the initial conditions were set. However, 
for gauge-independent tensor perturbations, a hypersur- 
face choice is irrelevant. Moreover, even for scalar per- 
turbations the entire integral can be evaluated in any 
of such traditional gauges as the synchronous, Newto- 
nian, comoving, or spatially flat gauges. Indeed, under a 

+ a'*, the integral trans- 



gauge transformation 5+ = 
forms as — 2njn A1 a^(r, x) + 2nin fJi a fJ i (Ti n , x- m ). The first 
of these two terms matches the transformation of the left 
hand side of eq. (A22). And the second one is negli- 
gible for a transformation of supcrhorizon perturbation 
from the uniform density gauge to any of the mentioned 
gauges [47]. 

In the integral of eq. (A22), the gradient of h^ v can be 
traded for time derivatives: 

n i h liV , i (r',x') = (-^-j - ±\ h^T^x'). 



The remaining integral d[iK{p, 1 x) can be calculated 
explicitly. However, even without its full calculation, it 
is easy to see that for x ^ — 1 



d^K{^,x) 



h( X ) [l + o(l)] 



(B3) 



Here is a proof of eq. (B2). This equation considers 
(1 - X 2 ) 2 /Ii d^iih{^)/(p - x) in the limit x -» -1 + 0. The 
most divergent at \ ~* ~ 1 contribution to f^dfi /i h{[i) / (/-t — x) 
arizes from the two discontinuous h(fi) terms (36). The products 
of either of the two corresponding integrals, 

dfi-^— = ln( X + l)[l + o(l)], 



dfi 



1 M - X 
Mln|l~M 2 | 
M - X 



1 



ln 2 ( X + l)[l + o(l)], 



and of the prefactor (1 
parently vanish as \ — > 



X 2 ) 2 in the considered expression ap- 
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Thus we can write eq. (31) as 

[{ X 2 -l)h!}' = -2R v h{ X ) [l + o(l)]. (B4) 
For x > — 1, the last equation is formally solved by 

h = -2R V [ X -P- r d X 2 h(X2) [1 + o(l)] + 

J-i xi - i J -i 

+ const. (B5) 

Given h of the general form (36), the above double inte- 
gral is continuous at x = — 1. Since the prefactor of h" in 
eq. (B4) vanishes at \ = ~ 1) the last constant in eq. (B5) 
may differ for x < — 1 an d X > — 1- The change of the 
constant provides a finite jump in h at x = — 1- This 
proves eq. (Bl). 

1. Subleading in 1/fcr term 

The result (B5) allows to establish a simple relation 
between the leading and subleading in 1/fcr terms that 



describe the subhorizon evolution of tensor modes during 
radiation domination. Evaluating the right derivatives of 
both sides of this equation at \ = — 1, wc find a relation 



h'(-l + 0) =R v h{-l + 0). (B6) 



Since h(\) is even, we have a similar equality at x = 1: 
h'(l - 0) = -R v h(l - 0). Remembering that h = for 
|x| > 1 and considering the Fourier components of h. we 
conclude that 



h^ d \kr) = A - Rv-p^r [i + o(i)] j . (B7) 



This relation applies to all orders of the R v expansion; 
its validity for 0{Ry) order is evident from eqs. (42) 
and (45). 
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